#include "Riostream.h"
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TMath.h"
#include "TRandom.h"
#include "TFormula.h"
#include "TSystem.h"
#include "TEnv.h"

#include "CommandLineInterface.hh"
#include "Barrel.hh"
#include "Annular.hh"
#include "Germanium.hh"
#include "SRIMloader.hh"

using namespace TMath;
using namespace std;

fstream Nfile;

ClassImp(Barrel);
ClassImp(Annular);
ClassImp(Germanium);

SRIM srim;

int main(int argc, char* argv[]){
	vector<char*> InputFiles;
	vector<char*> Histos;
	char* OutputFile = NULL;
	char* MBFile = NULL;
	char* cutfile = NULL;
	char* srimdir = NULL;
	char* SetFile = NULL;
	bool verbose = false;
	int tiltcorr = 0;
	int sim = 0;
	
	const Double_t pi = TMath::Pi();

	Int_t Np, Nd, Nt, Np_FB, Nd_FB, Nt_FB, Np_CD, Nd_CD, Nt_CD, tBe;
	Int_t par[3][5], ID[3][5];
	Int_t Apart, Atarget;
	Float_t Ein, thetain, phiin, xin, yin;
	Float_t Pxin, Pyin, Pzin;
	Float_t ttarget;
	Float_t theta, phi, thetadet, theta_cm[3][5];
	Float_t thetap[3][5], phip[3][5];
	Float_t theta_neu, phi_neu, dTheta, theta_neu2;
	Float_t x, y, z, r, alpha;
	Float_t xpad[4], ypad[4], zpad[4];
	Float_t Px, Py, Pz;
	Float_t vcm, vout[2], gamma[3][5], beta[3][5];
	Float_t Eex[3][5], Edet[3], Efor[3][5], EN, Q;
	Float_t dEcor, Ecor, Rafter;
	Double_t m[110], MBTheta[8][3][6], MBPhi[8][3][6];
	Double_t time[3][5], dt, dt2;
	Float_t CDThreshold[4];
	Float_t PT;
	Float_t EGamDS, EGamma, EGamma2;
	Int_t l, Ngamma;

	// Masses
	m[1] = 1.00782503207*931.494e3;
	m[2] = 2.01410177785*931.494e3;
	m[3] = 3.01604927767*931.494e3;
	m[4] = 4.00260325415*931.494e3;
	m[9] = 9.012182201*931.494e3;
	m[10] = 10.013533818*931.494e3;
	m[11] = 11.021657749*931.494e3;
	m[12] = 12.026920737*931.494e3; // 12Be
	m[20] = 19.999981315*931.494e3; // 20F
	m[21] = 20.999948951*931.494e3; // 21F
	m[22] = 22.002998815*931.494e3; // 22F
	m[23] = 22.994466904*931.494e3; // 23Ne
	m[107] = 106.905096820*931.494e3;

	CommandLineInterface* interface = new CommandLineInterface();

	interface->Add("-i", "inputfiles", &InputFiles);
	interface->Add("-o", "outputfile", &OutputFile);
	interface->Add("-M", "miniball angle", &MBFile);
	interface->Add("-C", "cut functions", &cutfile);
	interface->Add("-S", "SRIM files directory", &srimdir);
	interface->Add("-Set", "Setting file", &SetFile);
	interface->Add("-hi", "Histograms to be produced", &Histos);
	interface->Add("-t", "1 tiltcorr histos, 0 no tiltcorr", &tiltcorr);
	interface->Add("-v", "verbose", &verbose);
	interface->Add("-sim", "Set if you are analyzing a simulation",&sim);


	interface->CheckFlags(argc, argv);

	if(InputFiles.size() == 0 || OutputFile == NULL){
		cerr<<"You have to provide at least one input file and the output file!"<<endl;
		exit(1);
	}
	if(Histos.size() == 0){
		cerr << "No histograms will be filled!"<<endl;
		exit(1);
	}
	if(srimdir == NULL){
		cerr << "SRIM directory is needed!"<<endl;
		exit(1);
	}

	cout<<"input file(s):"<<endl;
	for(unsigned int i=0; i<InputFiles.size(); i++){
		cout<<InputFiles[i]<<endl;
	}

	cout<<"output file: "<<OutputFile<< endl;

	// Loading settings
	TEnv* Settings = new TEnv(SetFile);
	
	// Beam energy
	Ein = Settings->GetValue("E.beam",2.85e3);
	Ein = Ein*11;
	thetain = Settings->GetValue("theta.beam",0.);
	thetain = thetain/180.*pi;
	phiin = Settings->GetValue("phi.beam",0.);
	phiin = phiin/180.*pi;
	xin = Settings->GetValue("x.beam",0);
	yin = Settings->GetValue("y.beam",0.);

	ttarget = Settings->GetValue("t.target",10.);
	Atarget = Settings->GetValue("A.target",2);

	if (Atarget==12) m[12] = 12.000000000*931.494e3; // 12C

	// Pad settings
	for (Int_t ipad=0;ipad<4;ipad++){
		xpad[ipad] = Settings->GetValue(Form("x.pad.%d",ipad),0);
		ypad[ipad] = Settings->GetValue(Form("y.pad.%d",ipad),0);
		zpad[ipad] = Settings->GetValue(Form("z.pad.%d",ipad),0);
	}

	// A for particles stopped in dE
	Apart = Settings->GetValue("A.particle",2);

	// Generating SRIM graphs for energyloss calculations
	srim.Creategraphs(srimdir);
	Rafter = srim.RT[11]->Eval(Ein);
	Ein = srim.ET[11]->Eval(Rafter - ttarget/(2.*cos(thetain)));
	
	Pxin = sqrt(2*m[11]*Ein)*sin(thetain)*cos(phiin);
	Pyin = sqrt(2*m[11]*Ein)*sin(thetain)*sin(phiin);
	Pzin = sqrt(2*m[11]*Ein)*cos(thetain);
	
	TEnv *Angels = new TEnv(MBFile);
	for (int i=0; i<8; i++){
		for (int j=0; j<3; j++){
			for (int k=0; k<6; k++){
				MBTheta[i][j][k] = Angels->GetValue(Form("Theta.%d.%d.%d",i,j,k),0.0);
				MBPhi[i][j][k] = Angels->GetValue(Form("Phi.%d.%d.%d",i,j,k),0.0);
			}
		}
	}

	TChain* tr;
	tr = new TChain("caltr");
	for(unsigned int i=0; i<InputFiles.size(); i++){
		tr->Add(InputFiles[i]);
	}

	if(tr == NULL){
		cout << "could not find tree caltr in file " << endl;
		for(unsigned int i=0; i<InputFiles.size(); i++){
			cout<<InputFiles[i]<<endl;
		}
		return 3;
	}

	vector<Barrel> *FBarrel = new vector<Barrel>;
	vector<Barrel> *BBarrel = new vector<Barrel>;
	vector<Annular> *FCD = new vector<Annular>;
	vector<Germanium> *Miniball = new vector<Germanium>;
	long long EbisTime;
	long long T1Time;
	long long SuperCycleTime;

	tr->SetBranchAddress("ForwardBarrel",&FBarrel);
	tr->SetBranchAddress("BackwardBarrel",&BBarrel);
	tr->SetBranchAddress("ForwardCD",&FCD);
	tr->SetBranchAddress("Miniball",&Miniball);
	tr->SetBranchAddress("EbisTime",&EbisTime);
	tr->SetBranchAddress("T1Time",&T1Time);
	tr->SetBranchAddress("SuperCycleTime",&SuperCycleTime);

	TFile* outfile = new TFile(OutputFile,"recreate");

	if(outfile->IsZombie()){
		return 4;
	}

	CDThreshold[0] = 500.;
	CDThreshold[1] = 1300.;
	CDThreshold[2] = 500.;
	CDThreshold[3] = 800.;

	// Generating cut functions
	TFormula* For_cut[4][16];
	TFormula* Back_cut[4];
	TFormula* CD_cut[6][4];
	TFormula* CD_PT_cut[6][4];
	TFormula* Be_cut[4];

	TEnv *Cutfunc = new TEnv(cutfile);

	// Forward Barrel
	for (Int_t istrip=0;istrip<16;istrip++){
		for (Int_t icut=0;icut<4;icut++) For_cut[icut][istrip] = new TFormula(Form("For_cut_%d_%d",icut,istrip),Cutfunc->GetValue(Form("For.%d.%d",icut,istrip),"0*x"));
	}
	// 11 Beryllium
	Be_cut[0] = new TFormula("Be_cut_0",Cutfunc->GetValue("For.11","0*x"));
	Be_cut[1] = new TFormula("Be_cut_1",Cutfunc->GetValue("For.12","0*x"));

	// Backward Barrel
	for (Int_t icut=0;icut<4;icut++) Back_cut[icut] = new TFormula(Form("Back_cut_%d",icut),Cutfunc->GetValue(Form("Back.%d",icut),"0*x"));

	// CD
	for (Int_t idet=0;idet<4;idet++){
		for (Int_t icut=0;icut<6;icut++){
			CD_cut[icut][idet] = new TFormula(Form("CD_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD.%d.%d",icut,idet),"0*x"));
			CD_PT_cut[icut][idet] = new TFormula(Form("CD_PT_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD_PT.%d.%d",icut,idet),"0*x"));
		}
	}
	// 11 Beryllium
	Be_cut[2] = new TFormula("Be_cut_2",Cutfunc->GetValue("CD.11","0*x"));
	Be_cut[3] = new TFormula("Be_cut_3",Cutfunc->GetValue("CD.12","0*x"));

	// HISTOGRAMS
	// Prepare for vectorial histograms
	// Particle independent
	TH2F* For_dEvsEcor_Strip[4][16];

	// Proton
	TH1F* Proton_Eex[4];
	TH1F* Proton_Eex_Strip[4][16];
	TH1F* Proton_Eex_Angle[20];
	TH1F* Proton_Eex_320_Angle[20];
	TH1F* Proton_Eex_320_Angle_Background[20];
	TH1F* Proton_Eex_144_Angle[20];
	TH1F* Proton_Eex_144_Angle_Background[20];
	TH1F* Proton_Eex_144_Time_Angle[20];
	TH1F* Proton_Eex_144_Time_Angle_Background[20];
	TH1F* Proton_Eex_511_Angle[20];
	TH1F* Proton_Eex_511_Angle_Background[20];
	TH1F* Proton_Eex_511_Time_Angle[20];
	TH1F* Proton_Eex_511_Time_Angle_Background[20];
	TH1F* Proton_Eex_2100_Angle[20];
	TH1F* Proton_Eex_2100_Angle_Background[20];
	TH1F* Proton_Eex_2700_Angle[20];
	
	// Deuteron
	TH1F* Deuteron_Eex[4];
	TH1F* Deuteron_Eex_Strip[4][16];
	TH1F* Deuteron_Eex_Angle[20];
	TH1F* Deuteron_Eex_320_Angle[20];
	TH1F* Deuteron_Eex_320_Angle_Background[20];
	
	// Triton
	TH1F* Triton_Eex[4];
	TH1F* Triton_Eex_Strip[4][16];
	TH1F* Triton_Eex_Angle[20];
	TH1F* Triton_Eex_3300_Angle[20];
	TH1F* Triton_Eex_3300_comp_Angle[20];
	TH1F* Triton_Eex_2900_Angle[20];
	TH1F* Triton_Eex_2600_Angle[20];
	TH1F* Triton_Eex_6000_Angle[20];
	TH1F* Triton_Eex_6000_comp_Angle[20];
	
	// No punch through
	TH1F* NPT_Eex[4];
	TH1F* NPT_Eex_Strip[4][16];
	TH1F* NPT_Eex_Angle[20];
	TH1F* NPT_Eex_320_Angle[20];
	TH1F* NPT_Eex_320_Angle_Background[20];
	TH1F* NPT_Eex_511_Angle[20];
	TH1F* NPT_Eex_511_Angle_Background[20];
	TH1F* NPT_Eex_511_Time_Angle[20];
	TH1F* NPT_Eex_511_Time_Angle_Background[20];
	TH1F* NPT_Eex_2100_Angle[20];
	TH1F* NPT_Eex_2100_Angle_Background[20];
	TH1F* NPT_Eex_2600_Angle[20];
	TH1F* NPT_Eex_2700_Angle[20];
	TH1F* NPT_Eex_2900_Angle[20];
	TH1F* NPT_Eex_3300_Angle[20];
	TH1F* NPT_Eex_3300_comp_Angle[20];
	TH1F* NPT_Eex_3300_Angle_Background[20];
	TH1F* NPT_Eex_6000_comp_Angle[20];
	TH1F* NPT_Eex_6000_Angle_Background[20];


	// Creating histograms
	
	// Proton
	TH2F* Proton_EtotvsLabangletotal = new TH2F("Proton_EtotvsLabangletotal","E vs labangle for protons",180,0.,90.,500,0.,25.);
	TH2F* Proton_EtotvsCmangletotal = new TH2F("Proton_EtotvsCmangletotal","E vs Cmangle for protons",180,0.,180.,500,0.,25.);
	TH1F* Proton_Eextotal = new TH1F("Proton_Eextotal","Excitation energy for 12Be",1000,-2.5,7.5);
	TH1F* Proton_EGamma = new TH1F("Proton_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",8000,0,8.);
	TH1F* Proton_EDS = new TH1F("Proton_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",8000,0,8.);
	TH2F* Proton_EexvsEGamma = new TH2F("Proton_EexvsEGamma","Excitation energy vs gamma energy",200,-5.,15.,200,0,8.);
	TH2F* Proton_EexvsEDS = new TH2F("Proton_EexvsEDS","Excitation energy vs doppler corrected gamma energy",1000,-2.5,7.5,800,0,8.);
	TH2F* Proton_dtvsEGamma = new TH2F("Proton_dtvsEGamma","time diffens vs gamma energy",100,-1250.,1250.,8000,0,8.);
	TH2F* Proton_dtvsEDS = new TH2F("Proton_dtvsEDS","time diffens vs doppler corrected gamma energy",100,-1250.,1250.,8000,0,8.);
	TH2F* Proton_dtvsEex_144 = new TH2F("Proton_dtvsEex_144","time diffens vs excitation energy",100,-1250.,1250.,250,-2.5,7.5);
	TH2F* Proton_dtvsEex_511 = new TH2F("Proton_dtvsEex_511","time diffens vs excitation energy",100,-1250.,1250.,250,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_320 = new TH1F("Proton_Eex_Gamma_320","Excitation energy gated on .3-.34 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_400Less = new TH1F("Proton_Eex_Gamma_400Less","Excitation energy gated on .1-.4 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_400Less_Background = new TH1F("Proton_Eex_Gamma_400Less_Background","Excitation energy gated on 0.4-.8 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_144 = new TH1F("Proton_Eex_Gamma_144","Excitation energy gated on 0.14-0.147 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_144_Background = new TH1F("Proton_Eex_Gamma_144_Background","Excitation energy gated on 0.14-0.147 MeV",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_144_Time = new TH1F("Proton_Eex_Gamma_144_Time","Excitation energy gated on 0.14-0.147 MeV gammas + dT",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_144_Time_Background = new TH1F("Proton_Eex_Gamma_144_Time_Background","Excitation energy gated on 0.14-0.147 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_511 = new TH1F("Proton_Eex_Gamma_511","Excitation energy gated on 0.49-.52 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_511_Background = new TH1F("Proton_Eex_Gamma_511_Background","Excitation energy gated on 0.49-.52 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_511_Time = new TH1F("Proton_Eex_Gamma_511_Time","Excitation energy gated on 0.49-.52 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_511_Time_Background = new TH1F("Proton_Eex_Gamma_511_Time_Background","Excitation energy gated on 0.49-.52 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_2100 = new TH1F("Proton_Eex_Gamma_2100","Excitation energy gated on 2.07-2.15 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_2100_Background = new TH1F("Proton_Eex_Gamma_2100_Background","Background excitation energy gated on 2.07-2.15 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_2700 = new TH1F("Proton_Eex_Gamma_2700","Excitation energy gated on 2.68-2.8 MeV gammas",1000,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_2700_Background = new TH1F("Proton_Eex_Gamma_2700_Background","Background excitation energy gated on 2.07-2.15 MeV gammas",1000,-2.5,7.5);
	
	// Deuteron
	TH2F* Deuteron_EtotvsLabangletotal = new TH2F("Deuteron_EtotvsLabangletotal","E vs labangle for deuterons",180,0.,90.,500,0.,25.);
	TH2F* Deuteron_EtotvsCmangletotal = new TH2F("Deuteron_EtotvsCmangletotal","E vs Cmangle for deuterons",180,0.,180.,500,0.,25.);
	TH1F* Deuteron_Eextotal = new TH1F("Deuteron_Eextotal","Excitation energy",1000,-2.5,7.5);
	TH1F* Deuteron_EGamma = new TH1F("Deuteron_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* Deuteron_EDS = new TH1F("Deuteron_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",8000,0,8.);
	TH2F* Deuteron_EexvsEGamma = new TH2F("Deuteron_EexvsEGamma","Excitation energy vs gamma energy",200,-5.,15.,200,0,8.);
	TH2F* Deuteron_EexvsEDS = new TH2F("Deuteron_EexvsEDS","Excitation energy vs doppler corrected gamma energy",1000,-2.5,7.5,800,0,8.);
	TH1F* Deuteron_Eex_Gamma_320 = new TH1F("Deuteron_Eex_Gamma_320","Excitation energy gated on .3-.34 MeV gammas",1000,-2.5,7.5);
	
	// Triton
	TH2F* Triton_EtotvsLabangletotal = new TH2F("Triton_EtotvsLabangletotal","E vs labangle for tritons",180,0.,90.,500,0.,25.);
	TH2F* Triton_EtotvsCmangletotal = new TH2F("Triton_EtotvsCmangletotal","E vs Cmangle for tritons",180,0.,180.,500,0.,25.);
	TH1F* Triton_Eextotal = new TH1F("Triton_Eextotal","Excitation energy",1000,-2.5,7.5);
	TH1F* Triton_EGamma = new TH1F("Triton_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* Triton_EDS = new TH1F("Triton_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH2F* Triton_dtvsEGamma = new TH2F("Triton_dtvsEGamma","time diffens vs gamma energy",100,-1250.,1250.,8000,0,8.);
	TH2F* Triton_dtvsEDS = new TH2F("Triton_dtvsEDS","time diffens vs doppler corrected gamma energy",100,-1250.,1250.,8000,0,8.);
	TH2F* Triton_EexvsEGamma = new TH2F("Triton_EexvsEGamma","Excitation energy vs gamma energy",200,-5.,15.,200,0,8.);
	TH2F* Triton_EexvsEDS = new TH2F("Triton_EexvsEDS","Excitation energy vs doppler corrected gamma energy",1000,-2.5,7.5,8000,0,8.);
	TH1F* Triton_Eex_Gamma_2600 = new TH1F("Triton_Eex_Gamma_2600","Excitation energy gated on 3.3-3.5 MeV gammas",1000,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_2900 = new TH1F("Triton_Eex_Gamma_2900","Excitation energy gated on 3.3-3.5 MeV gammas",1000,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_3300 = new TH1F("Triton_Eex_Gamma_3300","Excitation energy gated on 3.3-3.5 MeV gammas",1000,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_3300_comp = new TH1F("Triton_Eex_Gamma_3300_comp","Excitation energy gated on 2.2-3.5 MeV gammas",200,0.,8.);
	TH1F* Triton_Eex_Gamma_6000_comp = new TH1F("Triton_Eex_Gamma_6000_comp","Excitation energy gated on 4.6-6.3.2 MeV gammas",1000,-2.5,7.5);
	
	// No Punch through
	TH2F* NPT_EtotvsLabangletotal = new TH2F("NPT_EtotvsLabangletotal","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH2F* NPT_EtotvsCmangletotal = new TH2F("NPT_EtotvsCmangletotal","E vs CM angle in for dE detector",180,0,180,600,0,15.);
	TH1F* NPT_Eextotal = new TH1F("NPT_Eextotal","Excitation energy",1e3,-2.5,7.5);
	TH1F* NPT_EGamma = new TH1F("NPT_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",8000,0,8.);
	TH1F* NPT_EDS = new TH1F("NPT_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",8000,0,8.);
	TH2F* NPT_dtvsEGamma = new TH2F("NPT_dtvsEGamma","time diffens vs gamma energy",100,-1250.,1250.,8000,0,8.);
	TH2F* NPT_EGammavsEGamma = new TH2F("NPT_EGammavsEGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.,800,0.,8.);
	TH2F* NPT_EexvsEGamma = new TH2F("NPT_EexvsEGamma","Excitation energy vs gamma energy",200,-5.,15.,200,0,8.);
	TH2F* NPT_EexvsEDS = new TH2F("NPT_EexvsEDS","Excitation energy vs doppler corrected gamma energy",200,-5.,15.,200,0,8.);
	TH2F* NPT_dtvsEGamma_320 = new TH2F("NPT_dtvsEGamma_320","time diffens vs gamma energy",100,-1250.,1250.,100,0,1.);
	TH2F* NPT_dtvsEGamma_320_Background = new TH2F("NPT_dtvsEGamma_320_Background","time diffens vs gamma energy",100,-1250.,1250.,100,0,1.);
	TH2F* NPT_dtvsEex_511 = new TH2F("NPT_dtvsEex_511","time diffens vs excitation energy",100,-1250.,1250.,250,-2.5,7.5);
	TH1F* NPT_Eex_Gamma = new TH1F("NPT_Eex_Gamma","Excitation energy gated on gammas",1000,-2.5,7.5);
	TH1F* NPT_EGamma_320 = new TH1F("NPT_EGamma_320","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* NPT_Eex_neutron = new TH1F("NPT_Eex_neutron","High excitation energy gated on neutron",100,-2.,8.);
	TH2F* NPT_EtotvsLabangle_neutron = new TH2F("NPT_EtotvsLabangle_neutron","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH2F* NPT_EtotvsLabangle_320_neutron = new TH2F("NPT_EtotvsLabangle_320_neutron","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH2F* NPT_EtotvsLabangle_2neutron = new TH2F("NPT_EtotvsLabangle_2neutron","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_neutron_Background = new TH1F("NPT_Eex_neutron_Background","High excitation energy gated on neutron",100,-2.,8.);
	TH2F* NPT_EtotvsLabangle_neutron_Background = new TH2F("NPT_EtotvsLabangle_neutron_Background","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_2neutron = new TH1F("NPT_Eex_2neutron","High excitation energy gated on two neutrons",100,-2.,8.);
	TH1F* NPT_Eex_2neutron_Background = new TH1F("NPT_Eex_2neutron_Background","High excitation energy gated on two neutrons",100,-2.,8.);
	TH1F* NPT_Eex_Gamma_320 = new TH1F("NPT_Eex_Gamma_320","Excitation energy gated on 0.31-0.33 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_320 = new TH2F("NPT_EtotvsLabangle_320","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_320_Background = new TH1F("NPT_Eex_Gamma_320_Background","Excitation energy gated on 0.3-0.31 og 0.33-0.34 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_320_Background = new TH2F("NPT_EtotvsLabangle_320_Background","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_320_neutron = new TH1F("NPT_Eex_Gamma_320_neutron","Excitation energy gated on 0.31-0.33 MeV gammas and neutrons",100,-2.,8.);
	TH1F* NPT_Eex_Gamma_320_neutron_Background = new TH1F("NPT_Eex_Gamma_320_neutron_Background","Excitation energy gated on 0.3-0.31 og 0.33-0.34 MeV gammas and neutrons",100,-2.,8.);
	TH1F* NPT_Eex_Gamma_511 = new TH1F("NPT_Eex_Gamma_511","Excitation energy gated on 0.31-0.33 MeV gammas",1000,-2.5,7.5);
	TH1F* NPT_Eex_Gamma_511_Background = new TH1F("NPT_Eex_Gamma_511_Background","Excitation energy gated on 0.3-0.31 og 0.33-0.34 MeV gammas",1000,-2.5,7.5);
	TH1F* NPT_Eex_Gamma_511_Time = new TH1F("NPT_Eex_Gamma_511_Time","Excitation energy gated on 0.31-0.33 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_511_Time = new TH2F("NPT_EtotvsLabangle_511_Time","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_511_Time_Background = new TH1F("NPT_Eex_Gamma_511_Time_Background","Excitation energy gated on 0.3-0.31 og 0.33-0.34 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_511_Time_Background = new TH2F("NPT_EtotvsLabangle_511_Time_Background","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_2100 = new TH1F("NPT_Eex_Gamma_2100","Excitation energy gated on 2.06-2.16 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_2100 = new TH2F("NPT_EtotvsLabangle_2100","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_2100_Background = new TH1F("NPT_Eex_Gamma_2100_Background","Excitation energy gated on 2.16-2.26 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_2100_Background = new TH2F("NPT_EtotvsLabangle_2100_Background","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_2600 = new TH1F("NPT_Eex_Gamma_2600","Excitation energy gated on 2.52-2.68 MeV gammas",1000,-2.5,7.5);
	TH1F* NPT_Eex_Gamma_2700 = new TH1F("NPT_Eex_Gamma_2700","Excitation energy gated on 2.69-2.79 MeV gammas",1000,-2.5,7.5);
	TH1F* NPT_Eex_Gamma_2900 = new TH1F("NPT_Eex_Gamma_2900","Excitation energy gated on 2.8-3.0 MeV gammas",1000,-2.5,7.5);
	TH1F* NPT_Eex_Gamma_3300 = new TH1F("NPT_Eex_Gamma_3300","Excitation energy gated on 3.25-3.5 MeV gammas",1000,-2.5,7.5);
	TH1F* NPT_Eex_Gamma_3300_comp = new TH1F("NPT_Eex_Gamma_3300_comp","Excitation energy gated on 2.5-3.5 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_3300_comp = new TH2F("NPT_EtotvsLabangle_3300_comp","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_3300_Background = new TH1F("NPT_Eex_Gamma_3300_Background","Excitation energy gated on 3.5-3.6 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_3300_Background = new TH2F("NPT_EtotvsLabangle_3300_Background","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_6000_comp = new TH1F("NPT_Eex_Gamma_6000_comp","Excitation energy gated on 4.2-6.2 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_6000_comp = new TH2F("NPT_EtotvsLabangle_6000_comp","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_6000_Background = new TH1F("NPT_Eex_Gamma_6000_Background","Excitation energy gated on 6.2-6.6 MeV gammas",1000,-2.5,7.5);
	TH2F* NPT_EtotvsLabangle_6000_Background = new TH2F("NPT_EtotvsLabangle_6000_Background","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH2F* NPT_ThetavsTheta_parvsneu = new TH2F("NPT_ThetavsTheta_parvsneu","Theta particle vs theta neutron",90,0.,90.,90,0.,90.);
	TH2F* NPT_ThetavsTheta_neuvsneu = new TH2F("NPT_ThetavsTheta_neuvsneu","Theta particle vs theta neutron",90,0.,90.,90,0.,90.);
	TH2F* NPT_ClusterCrystal_neuvsneu = new TH2F("NPT_ClusterCrystal_neuvsneu","Difference in cluster and crystal",20,-10.,10.,20,-10.,10.);
	TH2F* NPT_PhivsEex_Gamma_320_neutron = new TH2F("NPT_PhivsEex_Gamma_320_neutron","Delta phi vs E*",180,-180.,180.,100,-2.,8.);
	TH2F* NPT_PhivsEex_2neutron = new TH2F("NPT_PhivsEex_2neutron","Delta phi vs E*",180,-180.,180.,100,-2.,8.);

	// Pad histograms
	for (Int_t i=0; i<4; i++){
		// Proton
		Proton_Eex[i] = new TH1F(Form("Proton_Eex_%d",i),Form("Excitation energy in pad %f",i),1000,-2.,5.);
		
		// Deuteron
		Deuteron_Eex[i] = new TH1F(Form("Deuteron_Eex_%d",i),Form("Excitation energy in pad %f",i),1000,-2.,5.);
		
		// Triton
		Triton_Eex[i] = new TH1F(Form("Triton_Eex_%d",i),Form("Excitation energy in pad %f",i),1000,-2.,8.);
		
		// No Punch through
		NPT_Eex[i] = new TH1F(Form("NPT_Eex_%d",i),Form("Excitation energy in pad %f",i),1000,-5.,15.);
	
		// Pad+Strip histograms
		for (Int_t j=0; j<16; j++){
			// Particle independent
			For_dEvsEcor_Strip[i][j] = new TH2F(Form("For_dEvsEcor_Strip_%d_%d",i,j),Form("dE vs E for det %d and strip %d",i,j),500,0,20,500,0,20);


			// Proton
			Proton_Eex_Strip[i][j] = new TH1F(Form("Proton_Eex_Strip_%d_%d",i,j),Form("Excitation energy in pad %f",i),1000,-2.,5.);
		
			// Deuteron
			Deuteron_Eex_Strip[i][j] = new TH1F(Form("Deuteron_Eex_Strip_%d_%d",i,j),Form("Excitation energy in pad %f",i),1000,-2.,5.);
		
			// Triton
			Triton_Eex_Strip[i][j] = new TH1F(Form("Triton_Eex_Strip_%d_%d",i,j),Form("Excitation energy in pad %f",i),1000,-2.,8.);
		
			// No Punch through
			NPT_Eex_Strip[i][j] = new TH1F(Form("NPT_Eex_Strip_%d_%d",i,j),Form("Excitation energy in pad %f",i),1000,-5.,15.);
		}
	}

	// Angle histograms
	for (Int_t i=0; i<20; i++){
			// Proton
			Proton_Eex_Angle[i] = new TH1F(Form("Proton_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_320_Angle[i] = new TH1F(Form("Proton_Eex_320_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-7.5,7.5);
			Proton_Eex_320_Angle_Background[i] = new TH1F(Form("Proton_Eex_320_Angle_Background_%d",i),Form("Background excitaion for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_144_Angle[i] = new TH1F(Form("Proton_Eex_144_Angle_%d",i),Form("Excitation energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_144_Angle_Background[i] = new TH1F(Form("Proton_Eex_144_Angle_Background_%d",i),Form("Excitation energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_144_Time_Angle[i]= new TH1F(Form("Proton_Eex_144_Time_Angle_%d",i),Form("Excitation energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_144_Time_Angle_Background[i] = new TH1F(Form("Proton_Eex_144_Time_Angle_Background_%d",i),Form("Excitation energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
//			Proton_Eex_511_Angle[i] = new TH1F(Form("Proton_Eex_511_Angle_%d",i),Form("Excitation energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_511_Angle[i] = new TH1F(Form("Proton_Eex_511_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
//			Proton_Eex_511_Angle_Background[i] = new TH1F(Form("Proton_Eex_511_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_511_Angle_Background[i] = new TH1F(Form("Proton_Eex_511_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_511_Time_Angle[i] = new TH1F(Form("Proton_Eex_511_Time_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			Proton_Eex_511_Time_Angle_Background[i] = new TH1F(Form("Proton_Eex_511_Time_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
//			Proton_Eex_2100_Angle[i] = new TH1F(Form("Proton_Eex_2100_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_2100_Angle[i] = new TH1F(Form("Proton_Eex_2100_Angle_%d",i),Form("Excitation for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.,8.);
			Proton_Eex_2100_Angle_Background[i] = new TH1F(Form("Proton_Eex_2100_Angle_Background_%d",i),Form("Background excitation for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.,8.);
//			Proton_Eex_2700_Angle[i] = new TH1F(Form("Proton_Eex_2700_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),150,1.,4.);
			Proton_Eex_2700_Angle[i] = new TH1F(Form("Proton_Eex_2700_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 50.,(i+1)*3. + 50.),500,-2.5,7.5);
			
			// Deuteron
			Deuteron_Eex_Angle[i] = new TH1F(Form("Deuteron_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),500,-2.5,7.5);
			Deuteron_Eex_320_Angle[i] = new TH1F(Form("Deuteron_Eex_320_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),500,-2.5,7.5);
			Deuteron_Eex_320_Angle_Background[i] = new TH1F(Form("Deuteron_Eex_320_Angle_Background_%d",i),Form("Background excitaion for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),500,-2.5,7.5);
			
			// Triton
			Triton_Eex_Angle[i] = new TH1F(Form("Triton_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),500,-1.5,8.5);
			Triton_Eex_3300_Angle[i] = new TH1F(Form("Triton_Eex_3300_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,4.,7.);
			Triton_Eex_3300_comp_Angle[i] = new TH1F(Form("Triton_Eex_3300_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,2.,7.);
			Triton_Eex_2900_Angle[i] = new TH1F(Form("Triton_Eex_2900_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,4.,7.);
			Triton_Eex_2600_Angle[i] = new TH1F(Form("Triton_Eex_2600_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,4.,7.);
			Triton_Eex_6000_Angle[i] = new TH1F(Form("Triton_Eex_6000_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,4.,7.);
			Triton_Eex_6000_comp_Angle[i] = new TH1F(Form("Triton_Eex_6000_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 60.,(i+1)*3. + 60.),150,4.,7.);
		
		// No Punch through
		NPT_Eex_Angle[i] = new TH1F(Form("NPT_Eex_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_320_Angle[i] = new TH1F(Form("NPT_Eex_320_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_320_Angle_Background[i] = new TH1F(Form("NPT_Eex_320_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_511_Angle[i] = new TH1F(Form("NPT_Eex_511_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_511_Angle_Background[i] = new TH1F(Form("NPT_Eex_511_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_511_Time_Angle[i] = new TH1F(Form("NPT_Eex_511_Time_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_511_Time_Angle_Background[i] = new TH1F(Form("NPT_Eex_511_Time_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_2100_Angle[i] = new TH1F(Form("NPT_Eex_2100_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_2100_Angle_Background[i] = new TH1F(Form("NPT_Eex_2100_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_2600_Angle[i] = new TH1F(Form("NPT_Eex_2600_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),150,3.,8.);
		NPT_Eex_2700_Angle[i] = new TH1F(Form("NPT_Eex_2700_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),150,3.,8.);
		NPT_Eex_2900_Angle[i] = new TH1F(Form("NPT_Eex_2900_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),150,3.,8.);
		NPT_Eex_3300_Angle[i] = new TH1F(Form("NPT_Eex_3300_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),150,3.,8.);
		NPT_Eex_3300_comp_Angle[i] = new TH1F(Form("NPT_Eex_3300_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_3300_Angle_Background[i] = new TH1F(Form("NPT_Eex_3300_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),500,-2.,8.);
		NPT_Eex_6000_comp_Angle[i] = new TH1F(Form("NPT_Eex_6000_comp_Angle_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),150,3.,8.);
		NPT_Eex_6000_Angle_Background[i] = new TH1F(Form("NPT_Eex_6000_Angle_Background_%d",i),Form("Excitaion energy for angles between %f and %f",i*3. + 80.,(i+1)*3. + 80.),150,3.,8.);
		
	}

	// END HISTOGRAMS

	bool partfound;
	Double_t nentries = tr->GetEntries();
	Int_t nbytes = 0;
	Int_t status;

	Nd = 0;
	Np = 0;
	Nt = 0;
	Nd_FB = 0;
	Np_FB = 0;
	Nt_FB = 0;
	Nd_CD = 0;
	Np_CD = 0;
	Nt_CD = 0;
	tBe = 1;

	for(int i=0; i<nentries;i++){
		for (int idet=0;idet<3;idet++){
			for (int ie=0;ie<5;ie++){
				par[idet][ie] = 0;
				time[idet][ie] = 0;
				Eex[idet][ie] = 0;
				gamma[idet][ie] = 0;
				beta[idet][ie] = 0;
			}
		}
		status = tr->GetEvent(i);
		if(status == -1){
			cerr<<"Error occured, couldn't read entry "<<i<<" from tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<endl;
			return 5;
		}
		else if(status == 0){
			cerr<<"Error occured, entry "<<i<<" in tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<" doesn't exist"<<endl;
			return 6;
		}
		nbytes += status;


		partfound =false;

		// Forward barrel
		for(unsigned int j=0; j<FBarrel->size(); j++){
			if((*FBarrel)[j].GetMult() ==0) continue;

			if( (*FBarrel)[j].GetID() < 0 || (*FBarrel)[j].GetID() > 3 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". forward barrel detector id is wrong: "<< (*FBarrel)[j].GetID()<<endl;
				continue;
			}

			if( (*FBarrel)[j].GetStripNr()[0] < 0 || (*FBarrel)[j].GetStripNr()[0] > 15 ){
				cerr<<"Error in entry "<<i<<": "<<(*FBarrel)[j].GetID()<<". forward barrel detector strip nr is wrong: "<<(*FBarrel)[j].GetStripNr()[0]<<endl;
				continue;
			}
			
			partfound=true;
			time[0][j] = (*FBarrel)[j].GetTime();
			Edet[0] = (*FBarrel)[j].GetEdet();
			
			// Cylindrical coordinates
			x = 50*((*FBarrel)[j].GetStripPos()[0]-.5) + xin*cos(pi/2.*(*FBarrel)[j].GetID()) + yin*sin(pi/2.*(*FBarrel)[j].GetID());
			y = 29. + xin*sin(pi/2.*(*FBarrel)[j].GetID()) - yin*cos(pi/2.*(*FBarrel)[j].GetID());
			z = 8. + 3.125*((*FBarrel)[j].GetStripNr()[0] + gRandom->Rndm()); 
			// Moving each detector
			for (Int_t ipad=0;ipad<4;ipad++){
				if ((*FBarrel)[j].GetID()==ipad){
					x += xpad[ipad];
					y += ypad[ipad];
					z += zpad[ipad];
				}
			}

			r = sqrt(pow(y,2) + pow(x,2));
			alpha = atan(x/y);
			// Spherical coordinates
			theta = atan(r/z);
			phi = alpha + pi/2.*((*FBarrel)[j].GetID()+1);
			if (phi< 0) phi += 2*pi;
			if (phi> 2*pi) phi -= 2*pi;
			// Detector angle
			thetadet = acos(sin(theta)*cos(alpha));
			
			
			// No punchthrough
			PT = srim.EDL[Apart]->Eval(140./cos(thetadet));
			if (((*FBarrel)[j].GetEdet() < 400.) && ((*FBarrel)[j].GetRear() > 500.) && ((*FBarrel)[j].GetRear()<PT)) {
				
				par[0][j] = 4;
				ID[0][j] = (*FBarrel)[j].GetID();
				thetap[0][j] = theta;
				phip[0][j] = phi;
				
				// ASSUMING ONE PARTICLE TYPE!!!!
				
				// Energy loss in carbon foil
				Rafter = srim.RMylar[Apart]->Eval((*FBarrel)[j].GetRear());
				Efor[0][j] = srim.EMylar[Apart]->Eval(Rafter + 11.57/cos(thetadet));
				if (Atarget==107){
					Rafter = srim.RT[20]->Eval(Efor[0][j]);
					Efor[0][j] = srim.ET[20]->Eval(Rafter + ttarget/(2.*cos(theta)));
				}
				else {
					Rafter = srim.RT[Apart]->Eval(Efor[0][j]);
					Efor[0][j] = srim.ET[Apart]->Eval(Rafter + ttarget/(2.*cos(theta)));
				}
				Px = sqrt(2*m[Apart]*Efor[0][j])*sin(theta)*cos(phi);
				Py = sqrt(2*m[Apart]*Efor[0][j])*sin(theta)*sin(phi);
				Pz = sqrt(2*m[Apart]*Efor[0][j])*cos(theta);
				Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
				EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-Apart]);
				
				Eex[0][j] = Ein + Q - Efor[0][j] - EN;
				// CM Angles
				vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
				vout[0] = sqrt(2*Efor[0][j]/m[Apart])*cos(theta);
				vout[1] = sqrt(2*Efor[0][j]/m[Apart])*sin(theta);
				theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
				if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;

				// Relativistic factors
				gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-Apart] - m[Apart]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-Apart])/m[11+Atarget-Apart];
				beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
				
				
				// Filling histograms
				NPT_Eex[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
				NPT_Eex_Strip[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill(Eex[0][j]/1.e3);
				
				if (((*FBarrel)[j].GetID()==2)||((*FBarrel)[j].GetID()==3)){
					NPT_EtotvsLabangletotal->Fill(theta/pi*180.,(*FBarrel)[j].GetRear()/1.e3);
					NPT_EtotvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,(*FBarrel)[j].GetRear()/1.e3);
					NPT_Eextotal->Fill(Eex[0][j]/1.e3);
					for (Int_t iAngle=0; iAngle<20.;iAngle++)
					{
						if ((theta_cm[0][j]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 80./180.*pi))){
							NPT_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
							if (sim){
								NPT_Eex_320_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								NPT_Eex_511_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								NPT_Eex_2100_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								NPT_Eex_2600_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								NPT_Eex_2700_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								NPT_Eex_2900_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								NPT_Eex_3300_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								NPT_Eex_6000_comp_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
							}
						}
					}
				}
			}

			// Particle identification
			if (((*FBarrel)[j].GetEdet() > 400.) && ((*FBarrel)[j].GetRear() > 500.)){
				Ecor = (*FBarrel)[j].GetEdet() + (1. - cos(thetadet))*(*FBarrel)[j].GetRear();
				dEcor = (*FBarrel)[j].GetRear()*cos(thetadet);
				
				For_dEvsEcor_Strip[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill(Ecor/1.e3,dEcor/1.e3);

				// Protons
				if ((dEcor > For_cut[0][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[1][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
					par[0][j] = 1;
					ID[0][j] = (*FBarrel)[j].GetID();
					Np++;
					Np_FB++;
					thetap[0][j] = theta;
					phip[0][j] = phi;
					
					// Energy loss in carbon foil
					Rafter = srim.RMylar[1]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
					Efor[0][j] = srim.EMylar[1]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[1]->Eval(Efor[0][j]);
					Efor[0][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
					Px = sqrt(2*m[1]*Efor[0][j])*sin(theta)*cos(phi);
					Py = sqrt(2*m[1]*Efor[0][j])*sin(theta)*sin(phi);
					Pz = sqrt(2*m[1]*Efor[0][j])*cos(theta);
					Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
					EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
					
					Eex[0][j] = Ein + Q - Efor[0][j] - EN;
					// CM Angles
					vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
					vout[0] = sqrt(2*Efor[0][j]/m[1])*cos(theta);
					vout[1] = sqrt(2*Efor[0][j]/m[1])*sin(theta);
					theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
					if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
					// Relativistic factors
					gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-1])/m[11+Atarget-1];
					//gamma[0][j] = (Ein + (m[22] + m[2] - m[23] - m[1]) - Eex[0][j] - Efor[0][j] + m[23])/m[23];
					beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
			
					// Filling histograms
					Proton_Eex[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
					Proton_Eex_Strip[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill(Eex[0][j]/1.e3);
					
					if ((*FBarrel)[j].GetID()!=0){
						Proton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[0][j]/1.e3);
						Proton_EtotvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,Efor[0][j]/1.e3);
						Proton_Eextotal->Fill(Eex[0][j]/1.e3);
						
						for (Int_t iAngle=0; iAngle<20.;iAngle++)
						{
							if ((theta_cm[0][j]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 50./180.*pi))){ 
								Proton_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								if (sim){
									Proton_Eex_320_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Proton_Eex_144_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Proton_Eex_511_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Proton_Eex_2100_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Proton_Eex_2700_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								}
							}
						}
					}
				}

				// Deuterons
				if ((dEcor > For_cut[1][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[2][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
					par[0][j] = 2;
					ID[0][j] = (*FBarrel)[j].GetID();
					Nd++;
					Nd_FB++;
				
					thetap[0][j] = theta;
					phip[0][j] = phi;
					// Energy loss in carbon foil
					Rafter = srim.RMylar[2]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
					Efor[0][j] = srim.EMylar[2]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[2]->Eval(Efor[0][j]);
					Efor[0][j] = srim.ET[2]->Eval(Rafter + ttarget/(2.*cos(theta)));
					Px = sqrt(2*m[2]*Efor[0][j])*sin(theta)*cos(phi);
					Py = sqrt(2*m[2]*Efor[0][j])*sin(theta)*sin(phi);
					Pz = sqrt(2*m[2]*Efor[0][j])*cos(theta);
					Q = m[11] + m[Atarget] - m[2] - m[11+Atarget-2];
					EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-2]);
					Eex[0][j] = Ein + Q - Efor[0][j] - EN;
					// CM Angles
					vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
					vout[0] = sqrt(2*Efor[0][j]/m[2])*cos(theta);
					vout[1] = sqrt(2*Efor[0][j]/m[2])*sin(theta);
					theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
					if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
					// Relativistic factors
					gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-2] - m[2]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-2])/m[11+Atarget-2];
					//gamma[0][j] = (Ein + (m[22] + m[2] - m[22] - m[2]) - Eex[0][j] - Efor[0][j] + m[22])/m[22];
					beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
				
					// Filling histograms
					Deuteron_Eex[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
					Deuteron_Eex_Strip[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill(Eex[0][j]/1.e3);
					
					if ((*FBarrel)[j].GetID()!=0){
						Deuteron_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[0][j]/1.e3);
						Deuteron_EtotvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,Efor[0][j]/1.e3);
						Deuteron_Eextotal->Fill(Eex[0][j]/1.e3);
	
						for (Int_t iAngle=0; iAngle<20.;iAngle++)
						{
							if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))){
								Deuteron_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								if (sim){
								Deuteron_Eex_320_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								}
							}
						}
					}
				}

				// Tritons
				if ((dEcor > For_cut[2][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor)) && (dEcor < For_cut[3][(*FBarrel)[j].GetStripNr()[0]]->Eval(Ecor))){
					par[0][j] = 3;
					ID[0][j] = (*FBarrel)[j].GetID();
					Nt++;
					Nt_FB++;
					thetap[0][j] = theta;
					phip[0][j] = phi;
					// Energy loss in carbon foil
					Rafter = srim.RMylar[3]->Eval((*FBarrel)[j].GetEdet() + (*FBarrel)[j].GetRear());
					Efor[0][j] = srim.EMylar[3]->Eval(Rafter + 11.57/cos(thetadet));
					Rafter = srim.RT[3]->Eval(Efor[0][j]);
					Efor[0][j] = srim.ET[3]->Eval(Rafter + ttarget/(2.*cos(theta)));
					Px = sqrt(2*m[3]*Efor[0][j])*sin(theta)*cos(phi);
					Py = sqrt(2*m[3]*Efor[0][j])*sin(theta)*sin(phi);
					Pz = sqrt(2*m[3]*Efor[0][j])*cos(theta);
					Q = m[11] + m[Atarget] - m[3] - m[11+Atarget-3];
					EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-3]);
					
					Eex[0][j] = Ein + Q - Efor[0][j] - EN;
					// CM Angles
					vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
					vout[0] = sqrt(2*Efor[0][j]/m[3])*cos(theta);
					vout[1] = sqrt(2*Efor[0][j]/m[3])*sin(theta);
					theta_cm[0][j] = atan(vout[1]/(vout[0]-vcm));
					if (theta_cm[0][j] < 0) theta_cm[0][j] += pi;
					// Relativistic factors
					gamma[0][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-3] - m[3]) - Eex[0][j] - Efor[0][j] + m[11+Atarget-3])/m[11+Atarget-3];
					beta[0][j] = sqrt(1. - 1./(gamma[0][j]*gamma[0][j]));
						
					
					// Filling histograms
					Triton_Eex[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
					Triton_Eex_Strip[(*FBarrel)[j].GetID()][(*FBarrel)[j].GetStripNr()[0]]->Fill(Eex[0][j]/1.e3);
					
					if ((*FBarrel)[j].GetID()!=0){
						Triton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[0][j]/1.e3);
						Triton_EtotvsCmangletotal->Fill(theta_cm[0][j]/pi*180.,Efor[0][j]/1.e3);
						Triton_Eextotal->Fill(Eex[0][j]/1.e3);
						Triton_Eex[(*FBarrel)[j].GetID()]->Fill(Eex[0][j]/1.e3);
	
						for (Int_t iAngle=0; iAngle<20.;iAngle++)
						{
							if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))){
								Triton_Eex_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								if (sim){
									Triton_Eex_3300_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Triton_Eex_3300_comp_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Triton_Eex_2600_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Triton_Eex_2900_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Triton_Eex_6000_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
									Triton_Eex_6000_comp_Angle[iAngle]->Fill(Eex[0][j]/1.e3);
								}
							}
						}
					}
				}
			}
		}
		// End forward barrel


		// Gammas
		for(unsigned int j=0; j<Miniball->size(); j++){
			if((*Miniball)[j].GetCluID()<0 || (*Miniball)[j].GetCluID()>7){
				cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Cluster id is wrong: "<<(*Miniball)[j].GetCluID()<<endl;
				continue;
			}

			// Filling histograms
			for(unsigned int k=0; k<(*Miniball)[j].GetCrystal().size();k++){
				if( (*Miniball)[j].GetCrystal()[k].GetCryID() <0 || (*Miniball)[j].GetCrystal()[k].GetCryID() >2){
					cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Crystal id is wrong: "<<(*Miniball)[j].GetCrystal()[k].GetCryID()<<endl;
					continue;
				}
				

				l=0;
				// Gammas in coincidence with particles
				for (Int_t idet=0;idet<1;idet++){
					for (Int_t ie=0;ie<5;ie++){
						
						// Protons
						if ((par[idet][ie]==1)&&(ID[idet][ie]!=0)){
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							Proton_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Proton_EDS->Fill(EGamDS/1.e3);
							Proton_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Proton_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
							
							Proton_dtvsEGamma->Fill(25.*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Proton_dtvsEDS->Fill(25.*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),EGamDS/1.e3);
							if ((EGamDS<340)&&(300<EGamDS)){
								Proton_Eex_Gamma_320->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Proton_Eex_320_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<380)&&(340<EGamDS)){
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Proton_Eex_320_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<400)&&(100<EGamDS)){
								Proton_Eex_Gamma_400Less->Fill(Eex[idet][ie]/1.e3);
							}
							if ((EGamDS<800)&&(400<EGamDS)){
								Proton_Eex_Gamma_400Less_Background->Fill(Eex[idet][ie]/1.e3);
							}
							if (((*Miniball)[j].GetCrystal()[k].GetCore()<154)&&(134<(*Miniball)[j].GetCrystal()[k].GetCore())){
							  Proton_Eex_Gamma_144->Fill(Eex[idet][ie]/1.e3);
							  Proton_dtvsEex_144->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),Eex[idet][ie]/1.e3);
							  for(Int_t iAngle=0; iAngle<20.; iAngle++)
							    {
							      if((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_144_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
							    }
							  if (25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie])>-150){
							    Proton_Eex_Gamma_144_Time->Fill(Eex[idet][ie]/1.e3);
							    for (Int_t iAngle=0; iAngle<20; iAngle++)
							      {
								if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_144_Time_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3); 
							      }
							  }
							}
							if ((((*Miniball)[j].GetCrystal()[k].GetCore()<134)&&(124<(*Miniball)[j].GetCrystal()[k].GetCore()))||(((*Miniball)[j].GetCrystal()[k].GetCore()<164)&&(154<(*Miniball)[j].GetCrystal()[k].GetCore())))
							  {
							    Proton_Eex_Gamma_144_Background->Fill(Eex[idet][ie]/1.e3);
							    for (Int_t iAngle=0; iAngle<20.; iAngle++)
							      {
								if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_144_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
							      }
							    if (25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie])>-150)
							      {
								Proton_Eex_Gamma_144_Time_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.; iAngle++)
								  {
								    if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_144_Time_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								  }
							      }
							  }
							if (((*Miniball)[j].GetCrystal()[k].GetCore()<520)&&(490<(*Miniball)[j].GetCrystal()[k].GetCore())){
							  Proton_Eex_Gamma_511->Fill(Eex[idet][ie]/1.e3);
							  Proton_dtvsEex_511->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),Eex[idet][ie]/1.e3);
							  for (Int_t iAngle=0; iAngle<20.;iAngle++)
							    {
							      if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_511_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
							    }
							  if (25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie])>-150){
							    Proton_Eex_Gamma_511_Time->Fill(Eex[idet][ie]/1.e3);
							    for (Int_t iAngle=0; iAngle<20.;iAngle++)
							      {
								if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_511_Time_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
							      }
							  }
							}
							if ((((*Miniball)[j].GetCrystal()[k].GetCore()<535)&&(520<(*Miniball)[j].GetCrystal()[k].GetCore()))||(((*Miniball)[j].GetCrystal()[k].GetCore()<490)&&(475<(*Miniball)[j].GetCrystal()[k].GetCore()))){
								Proton_Eex_Gamma_511_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_511_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
								if (25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie])>-150){
									Proton_Eex_Gamma_511_Time_Background->Fill(Eex[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++)
									{
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_511_Time_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
							}
							if ((EGamDS<2160)&&(2060<EGamDS)){
								Proton_Eex_Gamma_2100->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_2100_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2260)&&(2160<EGamDS)){
								Proton_Eex_Gamma_2100_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_2100_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2800)&&(2650<EGamDS)){
								Proton_Eex_Gamma_2700->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 50./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 50./180.*pi))) Proton_Eex_2700_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3000)&&(2800<EGamDS)){
								Proton_Eex_Gamma_2700_Background->Fill(Eex[idet][ie]/1.e3);
							}
						}
						// Deuterons
						else if ((par[idet][ie]==2)&&(ID[idet][ie]!=0)){
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							Deuteron_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Deuteron_EDS->Fill(EGamDS/1.e3);
							Deuteron_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Deuteron_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
							if ((EGamDS<340)&&(300<EGamDS)){
								Deuteron_Eex_Gamma_320->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Deuteron_Eex_320_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<380)&&(340<EGamDS)){
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Deuteron_Eex_320_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
						}
						// Tritons
						else if ((par[idet][ie]==3)&&(ID[idet][ie]!=0)){
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							Triton_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Triton_EDS->Fill(EGamDS/1.e3);
							Triton_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Triton_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
							Triton_dtvsEGamma->Fill(25.*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Triton_dtvsEDS->Fill(25.*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),EGamDS/1.e3);

							if ((EGamDS<3500)&&(3300<EGamDS)){
								Triton_Eex_Gamma_3300->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_3300_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3500)&&(2200<EGamDS)){
								Triton_Eex_Gamma_3300_comp->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_3300_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3000)&&(2800<EGamDS)){
								Triton_Eex_Gamma_2900->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_2900_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2750)&&(2500<EGamDS)){
								Triton_Eex_Gamma_2600->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_2600_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<6200)&&(5800<EGamDS)){
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_6000_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<6200)&&(4600<EGamDS)){
								Triton_Eex_Gamma_6000_comp->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 60./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 60./180.*pi))) Triton_Eex_6000_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
						}
						// No punch through
						else if ((par[idet][ie]==4)&&(ID[idet][ie]!=0)){
							EGamma = (*Miniball)[j].GetCrystal()[k].GetCore();
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							theta_neu = MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]];
							phi_neu = MBPhi[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]];
							Ngamma = Miniball->size()*(*Miniball)[j].GetCrystal().size();
							if (Ngamma>1){
								for(unsigned int j2=0; j2<Miniball->size(); j2++){
									for(unsigned int k2=0; k2<(*Miniball)[j].GetCrystal().size();k2++){
										if (k<k2 || j<j2){
											EGamma2 = (*Miniball)[j2].GetCrystal()[k2].GetCore();
											NPT_EGammavsEGamma->Fill(EGamma/1.e3,EGamma2/1.e3);
										}
									}
								}
							}

							dt = 25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]);
							if ((ID[idet][ie]==2) || (ID[idet][ie]==3)){
								NPT_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
								NPT_EDS->Fill(EGamDS/1.e3);
								NPT_dtvsEGamma->Fill(25.*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
								NPT_Eex_Gamma->Fill(Eex[idet][ie]/1.e3);
								NPT_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
								NPT_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
								

								// Gamma gates
								if ((EGamDS<330)&&(310<EGamDS)){
									NPT_Eex_Gamma_320->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_320->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_320_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
									
									// gamma+neutron
									if (Ngamma>1){
										for(unsigned int j2=0; j2<Miniball->size(); j2++){
											for(unsigned int k2=0; k2<(*Miniball)[j].GetCrystal().size();k2++){
												if (k<k2 || j<j2){
													EGamma2 = (*Miniball)[j2].GetCrystal()[k2].GetCore();
													NPT_EGamma_320->Fill(EGamma2/1000.);
													NPT_dtvsEGamma_320->Fill(25*((*Miniball)[j2].GetCrystal()[k2].GetTime()-time[idet][ie]),(*Miniball)[j2].GetCrystal()[k2].GetCore()/1.e3);
													dt2 = 25*((*Miniball)[j2].GetCrystal()[k2].GetTime()-time[idet][ie]);
													if ((EGamma2<610&&590<EGamma2)||(EGamma2<850&&830<EGamma2)||(EGamma2<1020&&1010<EGamma2)){
														NPT_Eex_Gamma_320_neutron->Fill(Eex[idet][ie]/1000.);
														NPT_EtotvsLabangle_320_neutron->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
														NPT_PhivsEex_Gamma_320_neutron->Fill((phip[idet][ie]-phi_neu)/pi*180.,Eex[idet][ie]/1.e3);
													}
													if ((EGamma2<560&&540<EGamma2)||(EGamma2<810&&790<EGamma2)||(EGamma2<1080&&1070<EGamma2)){
														NPT_Eex_Gamma_320_neutron_Background->Fill(Eex[idet][ie]/1000.);
													}
													
												}
											}
										}
									}
								}
								if (((EGamDS<340)&&(330<EGamDS))||((EGamDS<310)&&(300<EGamDS))){
									NPT_Eex_Gamma_320_Background->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_320_Background->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_320_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								
								if (((*Miniball)[j].GetCrystal()[k].GetCore()<520)&&(490<(*Miniball)[j].GetCrystal()[k].GetCore())){
									NPT_dtvsEex_511->Fill(25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie]),Eex[idet][ie]/1.e3);
									NPT_Eex_Gamma_511->Fill(Eex[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_511_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
									if (25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie])>-150){
										NPT_Eex_Gamma_511_Time->Fill(Eex[idet][ie]/1.e3);
										NPT_EtotvsLabangle_511_Time->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
										for (Int_t iAngle=0; iAngle<20.;iAngle++)
										{
											if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_511_Time_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
										}
									}
								}
								if ((((*Miniball)[j].GetCrystal()[k].GetCore()<535)&&(520<(*Miniball)[j].GetCrystal()[k].GetCore()))||(((*Miniball)[j].GetCrystal()[k].GetCore()<490)&&(475<(*Miniball)[j].GetCrystal()[k].GetCore()))){
									NPT_Eex_Gamma_511_Background->Fill(Eex[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_511_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
									if (25*((*Miniball)[j].GetCrystal()[k].GetTime()-time[idet][ie])>-150){
										NPT_Eex_Gamma_511_Time_Background->Fill(Eex[idet][ie]/1.e3);
										NPT_EtotvsLabangle_511_Time_Background->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
										for (Int_t iAngle=0; iAngle<20.;iAngle++)
										{
											if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_511_Time_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
										}
									}
								}
								if (((EGamDS<2160)&&(2060<EGamDS))&&(((*Miniball)[j].GetCrystal()[k].GetCore()<2100)||(2140<(*Miniball)[j].GetCrystal()[k].GetCore()))){
									NPT_Eex_Gamma_2100->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_2100->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2100_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if (((EGamDS<2260)&&(2160<EGamDS))&&(((*Miniball)[j].GetCrystal()[k].GetCore()<2100)||(2140<(*Miniball)[j].GetCrystal()[k].GetCore()))){
									NPT_Eex_Gamma_2100_Background->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_2100_Background->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2100_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<2680)&&(2520<EGamDS)){
									NPT_Eex_Gamma_2600->Fill(Eex[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2600_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<2800)&&(2650<EGamDS)){
									NPT_Eex_Gamma_2700->Fill(Eex[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2700_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<3000)&&(2800<EGamDS)){
									NPT_Eex_Gamma_2900->Fill(Eex[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2900_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<3500)&&(3250<EGamDS)){
									NPT_Eex_Gamma_3300->Fill(Eex[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<3500)&&(2500<EGamDS)){
									NPT_Eex_Gamma_3300_comp->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_3300_comp->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<3600)&&(3500<EGamDS)){
									NPT_Eex_Gamma_3300_Background->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_3300_Background->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<6200)&&(4600<EGamDS)){
									NPT_Eex_Gamma_6000_comp->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_6000_comp->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_6000_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								if ((EGamDS<6600)&&(6200<EGamDS)){
									NPT_Eex_Gamma_6000_Background->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_6000_Background->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
									for (Int_t iAngle=0; iAngle<20.;iAngle++){
										if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_6000_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
									}
								}
								
								// Neutron gates.
								if ((EGamma<720&&690<EGamma)||(EGamma<836&&830<EGamma)||(EGamma<1045&&1030<EGamma)){
									NPT_Eex_neutron->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_neutron->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);

									// Neutron-Neutron events
									if (Ngamma==2){
										for(unsigned int j2=0; j2<Miniball->size(); j2++){
											for(unsigned int k2=0; k2<(*Miniball)[j].GetCrystal().size();k2++){
												if (k<k2 || j<j2){
													EGamma2 = (*Miniball)[j2].GetCrystal()[k2].GetCore();
													theta_neu2 = MBTheta[(*Miniball)[j2].GetCluID()][(*Miniball)[j2].GetCrystal()[k2].GetCryID()][(*Miniball)[j2].GetCrystal()[k2].GetSegID()[l]];
													if ((EGamma2<610&&590<EGamma2)||(EGamma2<850&&830<EGamma2)||(EGamma2<1020&&1010<EGamma2)){
														NPT_Eex_2neutron->Fill(Eex[idet][ie]/1000.);
														NPT_EtotvsLabangle_2neutron->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
														NPT_ThetavsTheta_parvsneu->Fill(thetap[idet][ie]/pi*180.,theta_neu/pi*180.);
														NPT_ThetavsTheta_neuvsneu->Fill(theta_neu/pi*180.,theta_neu2/pi*180.);
														NPT_ClusterCrystal_neuvsneu->Fill((*Miniball)[j].GetCluID()-(*Miniball)[j2].GetCluID(),(*Miniball)[j].GetCrystal()[k].GetCryID()-(*Miniball)[j2].GetCrystal()[k2].GetCryID());
														NPT_PhivsEex_2neutron->Fill((phip[idet][ie]-phi_neu)/pi*180.,Eex[idet][ie]/1.e3);
													}
													if ((EGamma2<560&&540<EGamma2)||(EGamma2<810&&790<EGamma2)||(EGamma2<1080&&1070<EGamma2)){
														NPT_Eex_2neutron_Background->Fill(Eex[idet][ie]/1000.);
													}
												}
											}
										}
									}
								}
								if ((EGamma<690&&670<EGamma)||(EGamma<860&&850<EGamma)||(EGamma<1055&&1045<EGamma)){
									NPT_Eex_neutron_Background->Fill(Eex[idet][ie]/1.e3);
									NPT_EtotvsLabangle_neutron_Background->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
								}
							}
						}
					}
				}
			}
		}
		// End gammas

		if(i%1000 == 0){
			cout<<setw(5)<<setiosflags(ios::fixed)<<setprecision(1)<<(100.*i)/nentries<<" % done\r"<<flush;
		}
	}

	// Writing histograms

	// Single histograms
	// Proton
	Proton_Eextotal->Write("",TObject::kOverwrite);
	Proton_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	Proton_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	Proton_EGamma->Write("",TObject::kOverwrite);
	Proton_EDS->Write("",TObject::kOverwrite);
	Proton_EexvsEGamma->Write("",TObject::kOverwrite);
	Proton_EexvsEDS->Write("",TObject::kOverwrite);
	Proton_dtvsEGamma->Write("",TObject::kOverwrite);
	Proton_dtvsEDS->Write("",TObject::kOverwrite);
	Proton_dtvsEex_144->Write("",TObject::kOverwrite);
	Proton_dtvsEex_511->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_400Less->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_400Less_Background->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_144->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_144_Background->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_144_Time->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_144_Time_Background->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_511->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_511_Background->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_511_Time->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_511_Time_Background->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_2100->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_2100_Background->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_2700->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_2700_Background->Write("",TObject::kOverwrite);
	
	// Deuteron
	Deuteron_Eextotal->Write("",TObject::kOverwrite);
	Deuteron_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	Deuteron_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	Deuteron_EGamma->Write("",TObject::kOverwrite);
	Deuteron_EDS->Write("",TObject::kOverwrite);
	Deuteron_EexvsEGamma->Write("",TObject::kOverwrite);
	Deuteron_EexvsEDS->Write("",TObject::kOverwrite);
	Deuteron_Eex_Gamma_320->Write("",TObject::kOverwrite);

	// Triton
	Triton_Eextotal->Write("",TObject::kOverwrite);
	Triton_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	Triton_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	Triton_EGamma->Write("",TObject::kOverwrite);
	Triton_EDS->Write("",TObject::kOverwrite);
	Triton_EexvsEGamma->Write("",TObject::kOverwrite);
	Triton_EexvsEDS->Write("",TObject::kOverwrite);
	Triton_dtvsEGamma->Write("",TObject::kOverwrite);
	Triton_dtvsEDS->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_2600->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_2900->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_3300->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_3300_comp->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_6000_comp->Write("",TObject::kOverwrite);

	// No punch through
	NPT_Eextotal->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	NPT_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	NPT_EGamma->Write("",TObject::kOverwrite);
	NPT_EDS->Write("",TObject::kOverwrite);
	NPT_dtvsEGamma->Write("",TObject::kOverwrite);
	NPT_EGammavsEGamma->Write("",TObject::kOverwrite);
	NPT_EexvsEGamma->Write("",TObject::kOverwrite);
	NPT_EexvsEDS->Write("",TObject::kOverwrite);
	NPT_dtvsEGamma_320->Write("",TObject::kOverwrite);
	NPT_dtvsEGamma_320_Background->Write("",TObject::kOverwrite);
	NPT_dtvsEex_511->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma->Write("",TObject::kOverwrite);
	NPT_EGamma_320->Write("",TObject::kOverwrite);
	NPT_Eex_neutron->Write("",TObject::kOverwrite);
	NPT_Eex_neutron_Background->Write("",TObject::kOverwrite);
	NPT_Eex_2neutron->Write("",TObject::kOverwrite);
	NPT_Eex_2neutron_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_320->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320_Background->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_320_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320_neutron->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320_neutron_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_511->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_511_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_511_Time->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_511_Time->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_511_Time_Background->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_511_Time_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_2100->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_2100->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_2100_Background->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_2100_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_2600->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_2700->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_2900->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_3300->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_3300_comp->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_3300_comp->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_3300_Background->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_3300_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_6000_comp->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_6000_comp->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_6000_Background->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_6000_Background->Write("",TObject::kOverwrite);
	NPT_ThetavsTheta_parvsneu->Write("",TObject::kOverwrite);
	NPT_ThetavsTheta_neuvsneu->Write("",TObject::kOverwrite);
	NPT_ClusterCrystal_neuvsneu->Write("",TObject::kOverwrite);
	NPT_PhivsEex_Gamma_320_neutron->Write("",TObject::kOverwrite);
	NPT_PhivsEex_2neutron->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_neutron->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_320_neutron->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_2neutron->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_neutron_Background->Write("",TObject::kOverwrite);

	// Pad dependent
	for (Int_t i=0; i<4; i++){
		//Proton
		Proton_Eex[i]->Write("",TObject::kOverwrite);
		
		//Deuteron
		Deuteron_Eex[i]->Write("",TObject::kOverwrite);
		
		//Triton
		Triton_Eex[i]->Write("",TObject::kOverwrite);
		
		// No punch through
		NPT_Eex[i]->Write("",TObject::kOverwrite);
	
		// Pad+strip dependent
		for (Int_t j=0; j<16; j++){
			//Particle independent
			For_dEvsEcor_Strip[i][j]->Write("",TObject::kOverwrite);
			
			//Proton
			Proton_Eex_Strip[i][j]->Write("",TObject::kOverwrite);
			
			//Deuteron
			Deuteron_Eex_Strip[i][j]->Write("",TObject::kOverwrite);
		
			//Triton
			Triton_Eex_Strip[i][j]->Write("",TObject::kOverwrite);
		
			// No punch through
			NPT_Eex_Strip[i][j]->Write("",TObject::kOverwrite);
		}
	}



	// Angle histograms
	for (Int_t i=0; i<20; i++){
		Proton_Eex_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_144_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_144_Angle_Background[i]->Write("",TObject::kOverwrite);
		Proton_Eex_144_Time_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_144_Time_Angle_Background[i]->Write("",TObject::kOverwrite);
		Proton_Eex_320_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_320_Angle_Background[i]->Write("",TObject::kOverwrite);
		Proton_Eex_511_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_511_Angle_Background[i]->Write("",TObject::kOverwrite);
		Proton_Eex_511_Time_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_511_Time_Angle_Background[i]->Write("",TObject::kOverwrite);
		Proton_Eex_2100_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_2100_Angle_Background[i]->Write("",TObject::kOverwrite);
		Proton_Eex_2700_Angle[i]->Write("",TObject::kOverwrite);
		Deuteron_Eex_Angle[i]->Write("",TObject::kOverwrite);
		Deuteron_Eex_320_Angle[i]->Write("",TObject::kOverwrite);
		Deuteron_Eex_320_Angle_Background[i]->Write("",TObject::kOverwrite);
		Triton_Eex_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_3300_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_3300_comp_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_2600_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_2900_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_6000_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_6000_comp_Angle[i]->Write("",TObject::kOverwrite);
		
		// No punch through
		NPT_Eex_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_320_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_320_Angle_Background[i]->Write("",TObject::kOverwrite);
		NPT_Eex_511_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_511_Angle_Background[i]->Write("",TObject::kOverwrite);
		NPT_Eex_511_Time_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_511_Time_Angle_Background[i]->Write("",TObject::kOverwrite);
		NPT_Eex_2100_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_2100_Angle_Background[i]->Write("",TObject::kOverwrite);
		NPT_Eex_2600_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_2700_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_2900_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_3300_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_3300_comp_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_3300_Angle_Background[i]->Write("",TObject::kOverwrite);
		NPT_Eex_6000_comp_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_6000_Angle_Background[i]->Write("",TObject::kOverwrite);

	}

	outfile->Close();
	delete tr;


	printf("Total: Np = %d, Nd = %d, Nt = %d, t = %f min \n",Np/Histos.size(),Nd/Histos.size(),Nt/Histos.size(),EbisTime*25.e-9/60.);
	printf("Pr. Minute: Np = %f, Nd = %f, Nt = %f \n",Np/Histos.size()/(EbisTime*25.e-9/60.),Nd/Histos.size()/(EbisTime*25.e-9/60.),Nt/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in FB: Np = %f, Nd = %f, Nt = %f \n",Np_FB/Histos.size()/(EbisTime*25.e-9/60.),Nd_FB/Histos.size()/(EbisTime*25.e-9/60.),Nt_FB/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in CD: Np = %f, Nd = %f, Nt = %f \n",Np_CD/Histos.size()/(EbisTime*25.e-9/60.),Nd_CD/Histos.size()/(EbisTime*25.e-9/60.),Nt_CD/Histos.size()/(EbisTime*25.e-9/60.));


	Nfile.open("Npar.txt", fstream::out | fstream::app);

	Nfile << Np/(EbisTime*25.e-9/60.) << "	" << Nd/(EbisTime*25.e-9/60.) << "	" << Nt/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Npar_FB.txt", fstream::out | fstream::app);

	Nfile << Np_FB/(EbisTime*25.e-9/60.) << "	" << Nd_FB/(EbisTime*25.e-9/60.) << "	" << Nt_FB/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Npar_CD.txt", fstream::out | fstream::app);

	Nfile << Np_CD/(EbisTime*25.e-9/60.) << "	" << Nd_CD/(EbisTime*25.e-9/60.) << "	" << Nt_CD/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Time.txt", fstream::out | fstream::app);

	Nfile << EbisTime*25.e-9 << endl;

	Nfile.close();
	
	return 0;

}

